JLS SAEM

data <- readRDS('data_sim_p_1000.rds')
for(x in names(data)) assign(x, data[[x]], envir = globalenv())
U_full <- data$U
rm(data)

S.data.time <- S.data$obs
S.data.time.log.sum <- sum(log(S.data.time))

model <- SAEM_model( 
  function(sigma2, ...) -n/(2*sigma2),
  function(phi1, phi2, phi3, ...) mean((Y - get_obs(phi1, phi2, phi3) )^2 ), 'sigma2',
  
  # === Variable Latente === #
  latent_vars = list(
    # === Non linear model === #
    latent_variable('phi', dim = G, size = 3, prior = list(mean = 'mu', variance = 'omega2'),
                    add_on = c('zeta(phi1 = phi1, phi2 = phi2, phi3 = phi3, ...)' )),
    
    # === S.data model === #
    latent_variable('b', prior = list(mean = 'barb', variance.hyper = 'sigma2_b'),
                    add_on = c('zeta(b = b, ...) +',
                               'sum(h$eval(b = b, ..., i = c(1,2)))' )),
    latent_variable('alpha', prior = list(mean = 'baralpha', variance.hyper = 'sigma2_alpha'),
                    add_on = c('zeta(alpha = alpha, ...) +',
                               'alpha*h$eval(alpha = alpha,..., i = 3)'))
  ),

  # === Paramètre de regression === #
  regression.parameter = list(
    regression_parameter('beta', 1, function(...) SPGD(5, theta0 = beta,
                                                      step = 0.05, lambda = 1/sqrt(N),
                                                      normalized.grad = T,
                                                      zeta.der.B, N, zeta.B, 
                                                      Z$alpha,  Z$phi1, Z$phi2, Z$phi3,Z$b) )
  )
)

# ---  Initialisation des paramètres --- #
parameter0 <- parameter %>% sapply(function(x) x* runif(1, 1.1,1.4))

parameter0$barb <- 10
parameter0$sigma2 <- 1
parameter0$baralpha <- 1
#===============================================#
load.SAEM(model)
oracle <- maximisation(1, do.call(S$eval, var.true), parameter, var.true)
#==============================================================================#

init.options <- list(x0 = list(phi = c(1,80,4), b = parameter0$barb, alpha = parameter0$baralpha), 
                     sd = list(phi = c(.05, 1.5, .5), b = 1, alpha = .1) )

SAEM.options <- list(niter = 200, sim.iter = 5, burnin = 190, 
                adptative.sd = 0.6)

saem

p <- 4
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta

SAEM.options$RDS <- paste0(params$rds_filename, '_', p)
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)

## [1] "SAEM execution time = 45min 44sec"

## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation

## [[1]]

## 
## [[2]]

## 
## [[3]]

Result of the SAEM-MCMC

	sigma2	mu.1	mu.2	mu.3	omega2.1	omega2.2	omega2.3	barb	baralpha
Real value	0.0025	0.9032	89.9575	5.0079	0.0039	35.9179	0.6948	30.0000	45.0000
Estimated value	0.0025	0.8979	90.0274	4.9032	0.0040	34.2305	0.5823	18.2806	30.2388
Rrmse	0.0101	0.0059	0.0008	0.0209	0.0098	0.0470	0.1619	0.3906	0.3280

p <- 20
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta

SAEM.options$RDS <- paste0(params$rds_filename, '_', p)
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)

## [1] "SAEM execution time = 46min 29sec"

## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation

## [[1]]

## 
## [[2]]

## 
## [[3]]

Result of the SAEM-MCMC

	sigma2	mu.1	mu.2	mu.3	omega2.1	omega2.2	omega2.3	barb	baralpha
Real value	0.0025	0.9032	89.9575	5.0079	0.0039	35.9179	0.6948	30.0000	45.0000
Estimated value	0.0025	0.9004	90.1263	4.9105	0.0044	34.6121	0.5214	17.7837	30.0382
Rrmse	0.0136	0.0031	0.0019	0.0194	0.1135	0.0364	0.2495	0.4072	0.3325

p <- 100
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta

SAEM.options$RDS <- paste0(params$rds_filename, '_', p)
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)

## [1] "SAEM execution time = 46min 12sec"

## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation

## [[1]]

## 
## [[2]]

## 
## [[3]]

Result of the SAEM-MCMC

	sigma2	mu.1	mu.2	mu.3	omega2.1	omega2.2	omega2.3	barb	baralpha
Real value	0.0025	0.9032	89.9575	5.0079	0.0039	35.9179	0.6948	30.0000	45.0000
Estimated value	0.0026	0.9025	90.1870	4.9241	0.0047	34.7273	0.5773	15.3291	27.0420
Rrmse	0.0189	0.0007	0.0026	0.0167	0.1872	0.0331	0.1690	0.4890	0.3991

p <- 1000
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta

SAEM.options$RDS <- paste0(params$rds_filename, '_', p)
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)

## [1] "SAEM execution time = 36min 26sec"

## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation

## [[1]]

## 
## [[2]]

## 
## [[3]]

Result of the SAEM-MCMC

	sigma2	mu.1	mu.2	mu.3	omega2.1	omega2.2	omega2.3	barb	baralpha
Real value	0.0025	0.9032	89.9575	5.0079	0.0039	35.9179	0.6948	30.0000	45.0000
Estimated value	0.0025	0.9037	90.0198	5.0017	0.0040	35.3809	0.6134	4.2174	8.2451
Rrmse	0.0046	0.0005	0.0007	0.0012	0.0277	0.0150	0.1172	0.8594	0.8168